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Multiple  State  Cox  Model.  toflglg  of  Survival 

The  rest  majority  of  reliability  analyses  for  components  and/or 
systeas  as suae  that  the  equipment  under  consideration  is  in  either  one  of 
two  states:  functioning  or  Failed.  In  aany  situations  one  can  aeasure  and 
distinguish  between  various  "levels  of  performance"  for  both  the  system 
and  its  components .  Current  literature  is  largely  devoted  to  the  proba¬ 
bilistic  aspects  of  both  static  and  dynamic  multistate  systeas.  Parallel 
Investigations  in  clinical  data  analysis  deal  with  reliability  under  the 
analogous  concept  of  survival  analysis. 

In  this  section  a  model  for  the  analysis  of  survival  data  of  a  staged 
failure  phenomena  with  competing  risks  is  considered.  The  model  is  a  mul¬ 
tiple  state  extension  of  the  proportional  hazard  rate  model  of  Cox  (1972). 

It  uses  the  theory  of  Markov  processes  as  applied  to  staged  failure  phenom¬ 
ena  described,  for  example,  in  Chiang  (1968).  Cox's  method  of  partial  like¬ 
lihood  vas  extended  to  Obtain  the  estimates  of  covarlate  effects.  Maximum 
likelihood  estimation  of  the  intensity  matrix  is  developed.  Using  the 
results  of  Aalen  (1975)>  the  consistency  of  the  estimator  for  the  integrated 
intensity  matrix  is  established.  The  asymptotic  variance  and  the  weak  con¬ 
vergence  to  normality  of  the  estimator  is  obtained.  An  illustration  apply¬ 
ing  the  model  to  a  multiple  transition  situation  is  given. 

1.  Introduction 

In  the  following,  we  consider  the  stochastic  modeling  of  staged 


failure  phenomena  using  the  theory  of  Harkov  processes.  Complicating 
the  statistical  picture  Is  a  set  of  observed  covariates  on  each  item 
sampled.  The  motivation  for  this  work  came  from  a  desire  to  solve 


certain  problems  that  arise  in  the  study  of  cancer  data.  Typically, 
for  each  patient  in  a  cancer  study,  there  exist  a  set  of  measured 
front  end  variables  (e.g.  initial  white  blood  cell  count,  age,  sex, 
etc.)  as  well  as  a  set  of  periodically  updated  medical  observations 
on  each  patient.  The  problem  is  to  properly  account  for  the  effect 
of  these  observed  covariates  on  the  transition  probabilities  between 
various  stages  of  the  disease.  This  situation  is  most  often  modeled 
using  just  two  states:  disease  and  death.  Cox  (1972)  considers 
this  situation  and  uses  the  parametric  model 

X  Jt)  •  \(t)  exp  ,  (1.1) 

where  \.(t)  is  the  failure  rate  (hazard  function)  for  the  ith  patient, 

If 

a.(t)  is  the  ith  patient's  observed  set  of  covariates,  and  is  the 

transpose  of  the  vector  representing  the  effect  of  the  covariates. 

The  function  X  (t)  is  a  reference  failure  rate  and  may  be  considered 
o 

as  the  failure  rate  of  a  patient  with  covariate  vector  x.(t)  *  o. 

If  covariates  are  measured  as  deviations  from  a  mean,  then  X^ft;  is 
an  "average"  failure  rate.  In  a  study  where  the  determination  of  the 
effect  of  covariates  is  of  paramount  importance,  then  the  method  of 
partial  likelihood  (Cox  1975)  is  a  very  useful  procedure  for  estima¬ 
ting  the  vector  6  in  that  it  results  in  a  formulation  independent  of 

the  unknown  X  (t). 

o 

In  this  paper,  we  consider  the  multiple  transition  extension  of 
a  parametric  covariate  model  in  an  analysis  of  life  data  (e.g.  the 
Cox  Model)  using  the  theory  of  Markov  processes  applied  to  staged 
failure  phenomena  as  described  for  example  in  Chiang  (1966).  The 


failed  states  (e.g.  death)  in  such  a  model  will  of  course  be  absorb' 
ing  states.  Competing  risk  models  can  be  handled  by  allowing  for 
more  than  one  death  state  where  each  death  state  corresponds  to  a 
different  cause  of  death.  For  example,  in  the  illustration  discuss¬ 
ed  in  Section  5,  patients  undergoing  porta-caval  shunt  surgery  to 
correct  hepatic  bleeding  can  at  any  particular  time,  either  reside 
in  one  of  two  treatment  states,  (alive-free  of  hepatitis,  or  alive- 
having  had  hepatitis)  or  in  one  of  three  absorbing  states  (death  from 
bleeding,  death  from  other  cause,  or  withdrawn- lost  to  follow-up). 

2.  STAGED  FAILURE  PHENOMENA:  MARKOV  PROCESSES 

We  consider  the  situation  where  each  patient  in  a  population 
must  exist  in  one  of  several  states.  The  states  will  be  comprised 
of  transient  states  (i.e.  illness  states  such  as  various  stages  of  a 
disease,  various  complications,  etc.)  and  absorbing  states  (i.e. 
death  states).  Complicating  the  analysis  of  many  clinical  studies 
is  the  censored  state,  which  here  may  be  treated  as  another  absorbing 
state.  The  Markov  assumption  is  that  the  probability  distribution 
over  future  states  a  patient  might  be  in,  is  completely  determined 
by  his  current  state.  Specifically  it  does  not  depend  upon  how  the 
patient  arrived  at  his  current  state  or  how  long  he  has  been  there. 
(This  statement  can  be  modified  if  the  state  space  is  expanded  to  in¬ 
clude  these  properties). 

Let  denote  the  set  of  mutually  exclusive  states  i,  iz  A(o) . 

In  our  applications  the  cardinality  of  A(o *  will  be  quite  small 
(five  for  the  liver  patient  example).  The  probability  of  being  in 


3 


state  j  at  tine  t  given  that  one  was  in  state  i  at  tine  t,t  <  t  is 


denoted  by  the  function  P.,(x,t).  The  fact  that  a  patient  must  occupy 
soae  state  for  each  time  t  and  that  the  states  are  mutually  exclusive 
is  characterized  by  the  relations 

7  P .  .(x ,t)  *  1  for  all  r <tj  and  for  all  i ,  (2.1a) 

. 

3 

and  P..(t,t)  *  6..,  (2.1b) 

I'd  •'3 

where  6 .  .  is  the  Kronecker  delta. 
id 

The  intensity  matrix  v^Jt)  is  given  by 

Vijlx)  ~  ~ Pis{x*t)  I  (2*2) 

t*T. 

From  (2.1a),  the  intensity  matrix  satisfies 

\  v..(t)  *  0  for  all  i,  t.  (2.3) 

•  * 


The  intensity  has  the  interpretation  of  being  proportional  to  the 
probability  of  transition  from  state  i  to  state  3  at  time  t. 
Specifically, 

P..  (t.t+b)  «  6..  +  v..  (t)&  +  o(b),  (2.4) 

%3  '  1 3  13 

where  o(L)  is  such  that 


11a  o(L)/L  *  0  o  r 

L+0  1  * 

Letting  P(x,t)  be  the  matrix  (P^^(xtt))  and  V(t)  be  the  matrix 

(v..(t))t  we  have  the  following  well  known  relationships  between  P 
x3 

and  V. 
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I.  Chapman- Kolmogrov  forward  differential  equation: 

3-  “  P(x,t)  V(t);  x<t,  (2.6) 


with  initial  condition  P(t,t)  »  I.  (2.7) 

II.  Chapman- Kolmogrov  backward  differential  equation: 


P(x,t)  -  -V(x)  P(x,t);  t <t 


with  "initial  condition" 


P(ttt)  «  J. 


(2.8) 


(2.9) 


The  products  PV  in  (2.6)  and  VP  in  (2.8)  are  the  usual  matrix  inner 
products.  The  integral  equation  form  of  (2.7)  to  (2.9)  is 

P(x,t)  -  I  +  J*  P  ( x  1 8 )  dB  (xt8),  (2.10) 


where 


B(x,t)  -  /*  V(a)  da. 
x 


(2.11) 


For  a  further  discussion  of  the  development  of  (2.6)  through  (2.9)  see 
Chiang  (1968),  pp.  116-120. 

In  the  usual  application  of  the  above,  the  intensity  matrix 

is  assumed  to  be  almost  everywhere  continuous  so  that  equa¬ 
tions  (2.6)  and  (2.8)  involve  no  irregularities.  It  will  be  seen 
later,  however,  that  useful  results  are  obtainable  if  we  let  v..(t) 


have  the  form 


n  .  . 
*3 


t«2 


ijl 


*3 


(2.12) 


where  6(  )  is  the  Dirac  delta  function  and  the  are  isolated 

points  on  the  real  line.  (In  particular,  this  is  the  natural  "non 
parametric  parameterization"  of  in  (1.1).  See  Cox  (1972)). 


S 


Equation  (2.12)  says  that  v*j  transitions  are  only  allowed  to  occur 


at  the  isolated  points  £-2, . . .  ,n  .  .. 

Computations  are  facilitated  by  ordering  the  for  all  iji. 


Let  the  ordered  be  relabeled  y.  with  y,<y  _<...<  y  . 

1  K  12  Ttl 


Let  P(y^,  y^)  be  the  probability  transition  matrix  for  the  times 
vk  m  vk  ~  c  w*th  o<e<yfc  -  It  can  be  shown  that 

Pfy”,y^;  is  the  identity  matrix  except  for  the  (i,j) th  off  diagonal 
element  and  the  ith  diagonal  element  which  are  given  by 


pn,vk-vk>  *  1  '-‘m1’ 


(2.12a) 


and 


Pii(]ik*vk} 


exp 


(2.12b) 


respectively.  The  Chapman -Kolmogrov  equations  imply  for. 


U*,~2  -*2  <vJc,  <vk,+l‘m-<vk0  -*2  <Vk0+l, 


k2 


that 


Pftjjtg)  *  1 1  p^k*^k^  ’ 

k^k  j 


(2.13) 


The  next  section  discusses  the  parametric  estimation  of  the  co¬ 
variate  effects.  The  following  section  discusses  the  nonparametric 
estimation  of  the  transition  probabilities  themselves  given  knowledge 
of  the  covariate  effects. 


3.  COVARIATE  PARAMETERIZATION  AND  LIKELIHOOD  ESTIMATION 


From  (2.7)  or  (2.8),  it  is  seen  that  the  transition  probability 
matrix  P  is  completely  determined  by  the  intensity  matrix  V.  The 
parameterization  of  the  covariate  effect  is  thus  facilitated  through 


a  parameterization  of  V.  Following  (1.1),  th'  Cox  model  parameteri¬ 
zation  for  the  pth  patient  for  a  transition  from  state  i  to  state  j 


is  taken  as 


(t)  exp  [at.  a  ft)],  i+j. 


(3.1) 


where  a.,  is  a  vector  of  length  Mt  the  number  of  covariates  in  the 
vector  Note  the  constraint  in  (3.1)  that  t  not  equal  j .  For 

the  case  the  parameterization  of  vi?.*  (t)  follows  from  (2.3). 


That  is, 


vtpJ  (t) 

ll 


- 1 


(3.2) 


Likelihood  estimation  in  both  the  case  in  which  the  are 

parameterized  and  the  case  in  which  they  are  left  arbitrary  is  des¬ 
cribed  in  Fertig,  Murthy,  et  al  (1979).  Here,  we  only  consider  the 
latter  situation. 

Let  (t)  be  the  right  continuous  state  function  for  the  pth 

patient.  We  observe  V(i)  on  [0,7.1  .  Here  we  take  7  to  be  a  non- 

P  P  P 

random  study  time  established  apriori  for  the  pth  patient. 

If  the  patient  leaves  the  study  because  of  death  or  is  censored 
unexpectedly,  he  will  be  observed  in  that  death  or  censored  state 
until  7p.  This  approach  explicitly  takes  into  account  a  censored 
state  that  may  be  entered  only  randomly.  Presumably,  the  transition 
probabilities  of  interest  are  those  that  would  exist  in  the  hypothe¬ 
tical  situation  of  no  censoring.  Estimation  of  these  probabilities 
is  the  thrust  of  the  Kaplan-Meier  estimator  in  the  three  state  problem 
(life-death-censored).  Fleming  (1978)  describes  a  multiple  state  ex¬ 
tension  to  the  K-M  estimator  that  can  be  used  in  the  current  situation. 
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The  state  functions  fc^fx)  are  each  seen  to  be  step  functions 
on  [0, Tp]  with  a  finite  number  of  jumps. 

Cox's  method  of  partial  likelihood  will  be  used  to  estimate  the 
in  (3.1).  In  this  situation,  one  may  specify  as  the  conditioning 
variates  (5 ^  in  Cox' s:  notation,  as  respectively  the  risk 

set  at  each  transition  time  (defined  below)  and  the  information  that 
a  specific  transition  occurred.  With  this  approach  we  need  to  record 
(1)  Ordered  event  times:  t^<...<tm.  (2)  Individuals  undergoing 

transition:  .  (3)  State  being  left:  ,x  .  (4)  State 

l  m  l  m 

being  entered:  y,,...,y  .  (5)  "Risk"  sets:  R  (t .R  (t  ).  The 

i  m  x  _  l  x m 

i  m 

sets  R  (t)  for  arbitrary  t  are  defined  by 
XJ 

Vt;  *  -*i>- 


Each  R  (t)  contains  as  members  all  patients  who  are  at  "risk”  of 

xo 

leaving  state  x .. 

0 

Cox's  partial  likelihood  is  then 


M 


L  = 


exp  (a'x  u  zi  ;/  E 
Xj  yj  Xo  ke.R 


exp  (a'x.  y.  zk}  ’  C3*3) 

x.(t  .)  9  *3 

0  3 


The  jth  factor  in  the  partial  likelihood  in  a  sense  measures  the 

likelihood  that  it  is  individual  i  .  who  undergoes  the  transition 

x  .+y  .  at  time  t.  given  that  some  individual  in  if  (t .)  must  undergo 
3  3  3  x  j  3 

this  transition.  (This  likelihood  is  written  for  the  case  in  which 
there  are  no  ties). 

A  useful  feature  of  (3.3)  is  that  it  can  be  factored  into  in¬ 
dependent  expressions  for  each  observed  transition  x+y.  Thus,  the 
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1 


estimation  of  o  for  given  (x,y)  can  be  performed  independently  of 

the  estimation  of  a  -  for  x+x'  or  y=y'.  This  factorization  also 

holds  for  the  information  matrix  based  on  (3.3).  Specifically  we 

have,  for  the  q th  component  of  a  , 

xy 


3  In  L 
3  a 


S-jL  m  l  C I  **  r<?;exp  (a'  «fcV  I  exp  (a'  z  )] 

xyq  t.cT  keF  (t.)  xy  k  kcF  (t.)  xy  k 


3  xy 


X3  ° 


x .  3 

3 


(3.4) 


where 


mh  i  ti  -  v }  • 


For  specified  x3y  the  expression  (3.4)  for  all  q  comprise  the  usual 
set  of  likelihood  equations  for  estimating  the  q  dimensional  parameter 
vector  of  the  Cox  model  when  only  two  states  x  and  y  (besides  the 
censoring  state)  are  possible. 

It  can  be  seen  from  (3.4)  that  all  off  diagonal  terms  of  the  in¬ 
formation  matrix  involving  different  transition  states  will  be  zero. 
Thus,  the  information  matrix  will  be  block  diagonal  with  the  (x,y) 


block  given  by 


T  P<? 

xy 


3^  In  L _  r  (  r 


ktR  (t.) 
X3  3 


zl^^zf^^  exp  (a'.,  z,) 
k  k  v  i3  k 


/  l  exp  (a.'. 
KzB  (t  J 
X3  J 


.  .  a. 

V3  k 


}  -  (  I 


ktHJt.l 

3 


k’.1  !a'xV  V  /  (J 


l  exP 

kcR(t.)  xy 

*  3  0 


z^P^  exp  (a'  a,)\ 

P  v  xy  k  \ 

2 

v  )j 


(3.5) 


Of  course,  we  only  need  consider  those  vectors  a  and  corres- 

xy 

ponding  J  ^  for  those  transitions  of  interest.  These  ideas  will  be 
xy 
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made  clearer  in  the  next  section. 

Estimation  of  the  transition  probabilities  require  an  estimation 
procedure  for  the  for  each  of  interest.  Kaplan  and 

Meier  (1958)  provide  the  nonparametric  estimate  in  the  three  state 
problem  with  no  covariates.  Fleming  (1978)  gives  the  procedure  for 
the  multi-state  case.  Cox  (1972)  provides  an  adaptation  of  the 
Kaplan-Meier  estimate  when  covariates  are  present.  The  estimator 
proposed  here  is  the  simple  extension  of  the  above  procedures  to  the 
case  where  covariates  exist  in  a  multi-state  environment.  It  is 
convenient  to  obtain  the  estimator  as  a  maximum  likelihood  estimator 
in  the  space  of  all  functions  whose  centers  form  a  set  with  no  cluster 
point.  This  procedure  provides  an  estimator  with  certain  intuitively 
nice  properties.  Specifically,  the  observed  conditional  frequencies 
of  the  sample  will  match  those  of  the  estimator's.  Moreover,  it  re¬ 
duces  to  the  Kaplan-Meier  estimate  when  there  are  only  three  states. 

It  is  recognized  that  the  fact  that  the  estimator  is  maximum  likeli¬ 
hood  may  provide  no  advantage  since  we  will  be  treating  a  case  where 
the  number  of  parameters  and  the  number  of  patients  on  study  grow 
at  the  same  rate.  Even  though  the  estimation  procedure  seems 
to  focus  on  the  intensity  matrix,  the  transition  probabilities  are  of 
primary  interest.  This  is  fortunate  since  the  presentation  of  a 
Dirac  delta  function  as  an  estimate  of  intensity  in  its  own  right 
would  be  quite  unsatisfactory  in  most  situations.  However,  the  transi 
tion  probabilities  essentially  involve  integrating  over  the  estimated 
intensities  and  therefore,  a  great  deal  of  smoothing  is  effected. 

This  corresponds  to  the  familiar  case  of  deriving  the  empirical 
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distribution  function  as  a  maximum  likelihood  estimator.  Its  den- 
sity  is  of  course  composed  of  delta  functions. 


With  the  above  considerations  in  mind  let  us  write  the  intensity 
for  the  pth  patient  as 


^  *« t  tah  V' 

for  i+j 


(3.6) 


where  A..  is  the  set  of  indices,  £,  for  which  a  transition  is  possible. 
**0 

The  set  .H.  A  .  .  which  forms  the  set  of  centers  of  the  Dirac  delta 
functions,  is  assumed  to  have  no  cluster  point.  For  a  given  clini¬ 
cal  trial  we  can  only  consider  the  £...  which  are  less  than  max  T  . 

xj  Jt  P  P 

Thus  we  may  take  .H.  A.,  as  a  finite  set. 

0  10 

The  data  that  are  observable  with  model  (3.6)  are  given  in  the 
following  table: 

Table  I:  Observable  Data  with  Model  (3.6) 


Event  Times:  ^ ijn 

Risk  Sets:  R .  (f ,  R.(£..  ) 
x  xjl  *  x  x  jn 

Transition  Sets:  S . .  ) 

xj  x o  *  t-j  xjn 


The  above  table  is  constructed  for  each  i  and  j  such  that  i+j.  The 
risk  set  R.(t)  is  the  set  of  all  patients  in  state  i  at  time  t~ . 

The  transition  set  S..(t)  is  the  set  of  all  patients  that  undergo  an 
x+j  transition  at  time  t.  Either  or  both  of  these  sets  may  be  empty 
at  any  t.  Thus  if  no  patient  undergoes  a  transition  from  x+j  at 
t  «  say,  then  #[s^.  As  *n  (2.12b)  the  probability 

that  a  patient  in  survives  the  x+j  transition  may  be  written 
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(ijf'  * p >  m  “»  [-*«!  laij  V1'  l3'7) 


The  probability  of  a  transition  is  one  minus  this  quality.  With 
(3.7)  it  is  simple  to  write  the  likelihood  of  the  data  in  Table  1 


n .  . 

■t.  IT  hr 

1-7  I  pESi, 


n-«p 


T 


pe*iatjV  ‘  sid<Kij i; 


.  (p) , 

e«p  (—fa) 


f3.8) 


where 


a.(p'*  *  a...  exp  faC.  a 

ij£  tj£  p 


(3.9) 


We  will  proceed  with  estimation  of  a...  in  the  case  that  the 

1-J  * 

a.,  vectors  are  known.  The  partial  likelihood  estimates  of  these 
vectors  may  be  used  since  these  estimates  were  derived  independently 
of  the  a .... 

To  estimate  the  a..0,  we  must  maximize  (3.8).  If  F.(£..0)  is 

vj  X  If  X-  J  X 

empty,  then  the  corresponding  a.  ..  cannot  be  estimated.  By  conven- 

X 

tion  we  take  it  to  be  zero.  If  /?.(£...)  is  non-empty  but  S.. (£...) 

%  X  tj  If  J  X 

is,  we  see  from  (3.8)  that  the  maximum  occurs  for  a"...  *  0.  Thus, 

tj  x 

we  need  consider  only  those  times  E.  ..  for  which  events  actually 

X 

occurred.  In  order  to  emphasize  this  fact  we  relabel  the  E...  as 

?'*/  x 

and  consider  them  "event”  times.  Since  we  are  estimating  in¬ 
tensities  but  are  primarily  interested  in  transition  probabilities, 
it  is  clear  from  the  decoupled  nature  of  the  log  of  (3.8)  that  those 
transitions  that  are  not  going  to  be  included  in  the  final  state 
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space  (e.g.  transitions  into  a  censored  state)  nay  be  ignored. 

We  are  now  reduced  to  the  case  where  for  each  t..0,1f(S..(t..a))>l 

and  R.(t.  .J  is  non-empty.  Setting  the  derivative  of  the  log  of 
x  tj  * 

(3.8)  with  respect  to  a^.^  to  zero,  we  have 


•* l  V  '-uji  e,p  *p;; 

1  -  exp  «XP  (aij  ZP)* 


p^i1  Un'  - s 


*  i  ^  t  •  •  «  i 
XQ  l 


exp  (a'  .  z  )  for  all  and  A.  (3.10) 

xj  p 


In  the  case  of  ’’ties"  at  t..#,  equation  (3.10)  may  be  solved  itera- 
tively  for  a....  In  the  important  case  of  no  ties,  we  have  the  follow- 

1»  J  A( 

ing  easily  solved  relation  for  a . . ^ : 


exp  f-y,  .xp  M» 

■  3  -  exp  fa'.  I  _  :>/  l 


XQ  p 


ijl  pcRJt^^) 


exp  (a; .  *  (3.ii) 


p 


where  is  that  patient  in  Sxo(txdl)' 

It  can  be  seen  from  (3.11)  that  in  the  case  of  no  covariate 


effects  »  0J,we  have 


•xp  ' 1 '  ~TiJ7 


(3.12) 


whore  is  the  number  of  patients  in  RJt^^).  Expression  (3.12) 

leads  to  the  usual  Kaplan-Meier  estimate  of  the  survival  function 
in  the  three-state  no  covariate  case.  The  estimator  given  in  (3.11) 
is  identical  to  that  proposed  by  Kalbfleisch  and  Prentice  (1973)  in 


the  three  state  problem  (life-death-censored).  They  develop  their 
estimator  by  considering  the  grouping  of  continuous  time  to  failure 
data. 


From  (3.11),  we  see  that  a^.^  is  a  function  of  the  covariates 
of  the  patients  in  the  risk  set  As  mentioned  earlier,  the 

estimates  will  in  practice  be  found  using  the  partial  likelihood 

A 

estimates  o,.t  of  the  a.,  vectors.  Having  estimated  the  a .  . 's  and 
the  we  may  write  an  estimate  of  the  intensity  matrix  for  a 

hypothetical  person  with  covariate  vector  a  as 


v . , (tiz) 


(o) 


iO 


A  ^  . 

ft;  a)  exp  fa'.  z)t  ifj 


v  . .  ( t;z)  «■  -  l  v  .  .(t;z) 

XZ  -  i  - 


(o; 


Because  each  v  .  .  is  a  sum  of  delta  functions,  we  have 

k2 


P(xtt;z)  m  |”[  Pf v“ 


(3.13) 


k-k. 


where  as  in  section  2  the  are  the  ordering  of  the  such  that 


Pj,  7  <  T  <  y.  <y,  +2<...<vi.  <  t  <  v*. 

Kj-i  -  *2~  K2+l 

A 

and  for  *s  ^e  identity  matrix  except  for  the 

ith  row  which  is  given  by 


*VV  V*; 


exp  **P  *))» 

A  a 

3 -exp  «*P  (al^. 


otherwise 


(3.14) 


These  ideas  are  implemented  in  section  5  with  an  example.  The 
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consistency,  asymptotic  variance  and  weak  convergence  of  the  in 
tegrated  intensity  matrix  estimator 

L*x)  "•/*  dt ' 

are  given  in  the  next  section. 


4.  CONSISTENCY  OF  THE  INTEGRATED  INTENSITY  MATRIX  ESTIMATOR 


The  notation  used  in  this  section  is  defined  in  Appendix  A. 
The  (xtj)  element  of  the  integrated  intensity  matrix  is  defined  by 


B .  .(ttb,z)  *  .  ,(t' ,z)  dt' 

xd  J  n 


xo 


(4.1) 


From  Section  2  the  maximum  likelihood  estimator  of  (4.1)  is 
given  by 


-  J*  W 


X3 


(t' Mz)h . .  dt' 
xj 


L  ^  ^  h.. 

/  h . .  \  JT7?~ 

y  ln  /j.  _ HL _  1  U* 


for  n .t .  .  >2  , 

X  XQV 

for  Vijr  s  2 


(4.2) 


Aalen  (1975)  presents  an  estimator  for  B..  of  the  form 

xj 


B.A.(ttLtz)  -  \z  v[OJ.  A(t',z)  h..  dt' 

h.. 

-  i  - 

{  iJr 


•t  C'Jo)  A 


(4.3) 


Noting  that  for  small  x,  it  can  be  intuitively  seen 

that  the  two  estimators  given  by  (4.2)  and  (4.3)  respectively  are 
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nearly  the  same  for  large  n...  In  Appendix  B,  this  point  is  used 
to  prove  the  following: 

Theorem  4.1 

a  p^u 

The  maximum  likelihood  estimate  B  (t,L,z)  is  consistent  for 
estimating  B(t,b,z).  i.e.. 


Plim  B  F~M(t,L,z)  -  B( t,L,z)  «  0 . 

N-”*  1  1 


*  F  M 

B  .  In  particular  we  have 
Theorem  4.2 


7 


t  A 


Thus,  the  asymptotic  properties  that  hold  for  B  also  hold  for 


/  g*. :(s)  dU..(a) 


IQ 


where  W..(z)  is  the  Weiner  process  and  /  g.-<8)  dw..(a)  is  a  normal 
process  with  independent  increments  and 

_ 'mUL _ 

£  /  g/(z)  P .* A0,8,z)  exp  (a*..  z)  dz 


j'H  8 


ij 


Theorem  4.3 


11m  Vat 
JIT-** 


/iT 


(V* 


-  / 


( t,L,z)  -  B^j(t,t, 


")] 


A  J  *rr*;  p... (o.a.z)  ||  • 

i  z  1  1 


•xp  fa'  .  z)  dz 
1*0 


5.  ILLUSTRATION  OF  MULTIPLE  TRANSITION  ESTIMATION 


In  order  to  illustrate  the  use  of  these  techniques  in  the  analysis 
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of  data  from  an  actual  clinical  trial,  we  have  drawn  data  from  a 
comparison  of  two  protocols  for  treatment  of  hepatic  hemorrhage. 

The  sample  consists  of  89  patients  who  were  admitted  to  a  local 
medical  center  with  hepatic  hemorrhaging  which  was  in  most  cases  a 
complication  of  alcholic  liver  disease.  The  patients  were  randomized 
into  one  of  two  treatment  protocols.  Patients  in  Group  1  received 
the  standard  medical  treatment  for  this  disorder,  which  consists  of 
whole  blood  transfusion  and  antibiotics.  Patients  in  Group  2  re¬ 
ceived,  in  addition  to  the  standard  treatment,  a  porta-caval  shunt 
operation  to  route  blood  around  the  hemorrhaging  portion  of  the  liver. 
There  was  a  lag  time  between  admission  to  the  study  and  randomization 
into  treatment  groups  which  ranged  from  0  to  22  months.  At  the  time 
of  randomization,  a  complete  blood  chemistry  workup  as  well  as  a 
disease  history  was  taken.  Demographic  data  for  this  clinical  trial 
are  contained  in  Fertig,  Murthy,  et  al  (1979). 

During  the  follow-up  period  occurrence  times  of  several  events 
were  noted.  The  relevant  events  for  the  purpose  of  this  example  are 
death  from  hepatic  hemorrhage,  death  from  other  cause,  withdrawal  re¬ 
sulting  in  loss  to  follow-up,  and  a  rise  in  blood  bilirubin  level 
over  3mg.  which  was  taken  as  an  indicator  of  hepatitis  onset.  These 
considerations  define  five  mutually  exclusive  states: 


State  1: 
State  2: 
State  3: 
State  4: 
State  5: 


Patient  alive,  no  hepatitis  since  randomization. 
Patient  alive,  hepatitis  infection  has  occurred. 
Dead  from  hepatic  hemorrhage. 

Dead  from  other  causes. 

Withdrawing  lost  to  follow-up. 
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Since  State  1  can  never  be  reentered,  and  since  states  3  through 
S  are  absorbing  states,  the  only  allowable  transitions  are  1-2,  1-3, 
1-4,  1-5,  2-3,  2-4,  and  2-5.  All  patients  entered  the  study  in  state 
1. 

The  above  classification  of  states  represents  a  competing  risk 
model  when  there  is  the  risk  of  dying  from  hepatic  hemorrhaging  and 
the  risk  of  dying  from  other  causes.  The  primary  interest  in  this 
study  is  in  the  effect  of  treatment  on^the  probability  of  dying  from 
hepatic  hemorrhaging.  Since  the  treatment  for  each  of  the  patients 
in  this  study  involved  whole  blood  transfusions,  there  was  a  large 
incidence  of  hepatitis,  the  presence  of  which  could  alter  the  effec¬ 
tiveness  of  the  porta-caval  shunt  operation.  Thus,  it  is  important 
to  be  able  to  estimate  the  effect  of  treatment  on  the  1-3  transition 
intensity  (no  hepatitis  -*•  hepatic  hemorrhage)  separately  from  its 

effect  on  the  2-3  intensity  (hepatitis  -*•  hepatic  hemorrhage).  The 

0 

multiple  transition  model  presented  in  this  paper  allows  for  just 

this  separate  estimation  to  be  performed. 

To  this  end,  we  define  the  fiTSt  covariate,  z*2\  such  that 

3^*0  for  those  patients  undergoing  standard  treatment,  and  z^22  =  l 

for  those  patients  undergoing  a  porta-caval  shunt  operation.  For 

(2) 

illustrative  purposes,  we  consider  a  second  covariate,  z  ,  defined 
as  0  for  those  patients  whose  time  to  randomization  was  less  than 
the  median  time,  and  as  1  for  those  patients  whose  time  to  randomiza¬ 
tion  was  greater  than  the  median. 

The  total  number  of  transitions  for  all  69  patients  was  142.  A 
selected  subset  of  these  transitions,  together  with  covariate 
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information  is  presented  in  Table  1  below. 


TABLE  1 

Selected  Events  for  Porta-Caval  Shunt  Clinical  Trial 


Treatment  Leg  Time 

Event  From  To  Covariate  Covariate 

Event  Patient  Time  State  State  Value,,,  Value,.,  D  , . 

No.  I.D.  (Weeks)  i  J_ _ a  (2J  s  (Z) 


1  33 

2  69 

3  75 

4  69 


0.286  1 

0.571  1 

0.714  1 

2.000  2 


4  1 

2  0 

4  1 

3  0 


1  (1,2, ...,89} 

0  (33) 

0  (33,69) 

0  (69) 


70 


47  84.286  2  4 


0  (7,13,18,20, 

23,29,31,47, 
54,56,62,64  , 
70,72,76,80, 
81,85} 


139  12 

140  9 

141  1 

142  3 


476.143  1 

522.143  1 
523.428  1 
525.857  1 


5  1 

5  1 

5  1 

5  1 


1  (13,9,12) 

0  (1,3,9) 

0  (1,3) 

1  (3} 


J 


Transitions  i+St  i*l,  2  represent  loss  to  follow-up.  These  tran¬ 
sitions  supply  no  information  pertinent  to  the  study  of  treatment 
effect.  State  5  is  thus  not  treated  as  a  competing  risk.  Rather,  the 
time  that  a  patient  undergoes  in  t-*-5  transition  is  the  time  he  is 
deleted  from  P^(t)  (the  risk  set  of  individuals  who  may  undergo  a 
transition  from  state  i).  This  is  formally  done  by  setting 
for  all  at  which  these  transitions  have  occurred.  In  this  case 
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the  matrix  whose  elements  are  defined  by  equation  (3.14)  reduces  to 
the  identity,  and  the  matrix  of  transition  probabilities  remains  un¬ 
changed  at  time  t 

The  raw  estimates  of  the  parameters  a^.  were  obtained  by  the 
partial  likelihood  methods  discussed  in  Section  3  and  are  displayed 
in  Table  2  along  with  their  standard  scores,  obtained  by  dividing 
each  estimate  by  its  estimated  asymptotic  standard  error. 


TABLE  2  / 

Covariate 

Parameter  Estimates 

From 

To 

Treatment  Effect 

id 

Lag  Time 

Effect  a(.V 

i«7 

State 

State 

Estimate 

Est / S .E . 

Estimate 

Est/S .E . 

1 

2 

-0.0139 

-0.05 

0.0486 

0.17 

1 

3 

<•00 

- 

0.2432 

0.33 

1 

4 

0.0687 

0.11 

0.2837 

0.46 

2 

3 

-1.9088 

-2.43 

0.6218 

1.07 

2 

4 

0.8275 

1.65 

0.6937 

1.55 

*  (j  \ 

The  estimate  of  *  -«  resulted  because  all  the,  1-3  transi¬ 

tions  (eight  of  them)  occurred  among  people  in  the  "standard  treat¬ 
ment"  group.  The  permutation  test  of  the  hypothesis  that  the  eight 
transitions  are  evenly  distributed  between  the  treatment  and  no  treat¬ 
ment  groups  rejects  this  hypothesis  at  the  0.01  significance  level. 

Using  the  techniques  of  Section  3,  the  matrix  of  transition 
probabilities  was  estimated  at  each  event  time.  For  example,  consider 
the  first  event  in  Table  1.  A  patient  with  covariate  a  »  (1,1)  en- 

tered  state  4  from  state  1.  From  Table  2,  we  can  compute  a.'  s  * 

Pj42 

0.3623.  At  the  time  of  this  event,  all  the  patients  were  in 
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state  1,  so  all  were  at  risk.  The  value  of  £  exp  (a'24*^)  for  this 
risk  set  is  107.283.  The  resulting'estimate  of  Zj4i  *rom  ©Quati 
(3.11)  is  *241  *  "*n 


on 


0. 3S2S/107.  283)yg0. 3S23 


« 0,009383 .  The 


probability  transition  matrix,  P .  .(0,t;z)  is  the  identity  matrix  for 

t  j 

t<t (22*20.286.  At  t*20.  286,  it  is  updated  according  to  equations 
(3.13)  and  (3.14).  For  the  case  z*(0,04,  we  see  from  (3.14),  that 


0.9907 

0 

0 

0.00934 

0 

fa 

1 

0 

0 

0 

P(0,  t.;z)*P(t~,  t.;z*(0,0))  = 

0 

0 

1 

0 

0 

1  *  1 

0 

0 

0 

1 

0 

0 

0 

0 

0 

1 

“ 

A 

ft  *  P  . 

- 

ft  ^  ? 

where  the  P^4  element  is  given  by  1  -[exp  -0.009383  e  14  ]  which 

reduces  to  0.00934  for  z*(0,0). 


As  an  illustration  of  the  updated  procedure,  we  consider  the 
event  at  time  t (? q2*84 .  286 .  Without  specifying  the  details,  we  have 
*24t  -0-  01987.  From  this  we  can  compute  P^At~7Q,  t70;z)  from 
equation  (3.13).  The  non  identity  portion  of  this  matrix  is 

A  ^ 

P24(t70*  t70iz)  *  J_exP  [ -0.  01987  e  a24Z  ]  *  0.01967  for 
z=  ( 0,  0)  and  P22(t~70 3  t70*Z *  =  0'980SS- 


The  updating  of  P(0,t;z)  to  t*t?Q  is  done  by  noting 

P(0tt7QiZ)  *  P(0,tggJZ)  P  (  t  Q g,  t  7  q!  Z  )  P  (  t  7  Q,  t7QjZ). 


But  P(t^g,  t~0;z)  *  I.  Thus,  for  **(0,0), 
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0.403  0. 223  0.264  0.111 


where  P( 0 ,t  gg;z)  is  as  indicated  above. 

We  see  from  this,  for  example,  that  the  probability  a  patient 
with  covariate  value  z-(ot0)  remains  in  state  1  for  t7Q=84.Z  weeks 
is  estimated  as  0.403.  Note  that  the  fifth  state,  lost-to-follow-up, 
has  zero  entry  probability  from  states  1-4.  This  is  a  consequence 
of  setting  for  and  is  in  keeping  with  the  Kaplan  - 

Meier  treatment  of  the  censored  state. 

We  present  in  Figures  1  and  2  plots  of  ?13(0,t;z)  and  P  (o,t;z) , 
respectively  for  all  values  of  a.  A  literal  interpretation  of  Figure 
1  is  that  a  person  in  the  surgical  treatment  groups  is  much  less 
likely  to  die  from  hepatic  hemorrhage,  than  a  person  in  the  standard 
treatment  group,  which  is  the  expected  result.  This  large  difference 


22 


is  primarily  due  to  the  infinite  estimate  for  There  is 

essentially  no  difference  that  can  'be  attributed  to  lag  time  to 
randomization. 


From  Figure  2  we  see  a  person  in  the  surgical  group  is  much 
more  likely  to  eventually  die  from  other  causes  than  a  person  in  the 
non-treatment  group.  This  apparent  negative  consequence  of  treatment 
could  easily  be  accounted  for  by  the  fact  that  those  patients  who 
were  not  treated  had  a  higher  incidenc^  of  death  from  hepatic  hemorr¬ 
haging  and  thus  were  not  available  to  die  from  "other  causes".  In¬ 
spection  of  Table  2  indicated,  however,  a  very  marginal  positive 
a^.  Though  not  statistically  significant,  an  of  0.8275  is  an 

indication  that  treatment  enhances  the  hazard  of  dying  from  other 


FIGURE  t 


Probability  of  death  due  to  hepatic  heaorrbaee  in  the  preaeace  of 
all  other  causes  of  death  versus  tiae  oa  study. 


Elapsed  tiae  in  Study:  b'eeks 


FIGURE  2 

Probability  of  death  due  to  other  causes  in  the  presence  of  all  risks 
of  death  versus  tiae  on  study. 
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APPENDIX  B 

Proof  of  Consistency  and  Weak  Convergence 


We  will  now  prove  the  consistency  of  our  estimator  for  the  integrated 
intensity  matrix 

B..(t,L,»)  -  /*  v.  .  ( tj  z )  dt '  CD 

ij  A  tj 


Assuming  the  a.,  vectors  are  known,  our  estimate  for  (1)  denoted  by 
B  F~M (t,htz)  is  given  by 


B  F'M(t,L,z) 


F-M 
( t,z) 


h .  .  dt 


*  l 


In 


/  h  .  . 

*  ijr 


n  .  h  .  . , 

1ft  •  •  tr  J  ^  •  • 

ijr  tjr 


-  0 


for  n  ..  >2 

it .  . 
ijr 


for  n.,  <2 

ltijr 


(2) 


We  will  prove  the  consistency  of  our  estimator  (2)  for  estimating  (1) 
by  proving  that  our  estimator  is  equivalent  in  probability  to  Aalen's 
consistent  estimator  (1975)  for  estimating  (1).  Now  Aalen's  estimator 
for  estimating  (1)  is  given  by 


BA(ttb,z)  =  J*  v °..A(t'z)  h..  dt' 

*  *»  n  9  n  n 


u 


tj 

h  .  • 

.  .bL 


L<t .  .  <t  n  . ,  h .  .. 

•  tjr-  it .  .  ijt .  . 

ijr  *  ijr 


(3) 


We  will  now  demonstrate  that  the  difference  between  (2)  and  (3)  will 
converge  to  zero  in  probability.  To  simplify  the  proof  we  assume 
the  following: 
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Inf 

it  j, xeZ 


exp 


>  0 


A 


2  •* 
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Sup 

i,jt  *cZ 


exp 


<  ® 


27 


A3  *’  iil  Finit>2'  *  •  •  *nat>2*  v0<^tST)  -  2 
A.  •*  ;S^«7  (t)  exp-'  (a'.  .  z) „  -  K  <  <*  . 

9  IjJZtZ  IJ  tj 

The  proposition  to  be  established  is  then  given  by  the  following 
theorem  and  its  proof. 


Theorem: 


*  A(t,L,z)  -  B  F~M(ttL,z)  -  0 


for  all  t  such  that 


L<t<T. 


Proof : 


For  given  e>0,  8>0 ,  in  view  of  choose  N  (T)  such  that  for  N>NQ(t) 
P  (n  . .  >2,  . . .  ,n  .  >2  y.  3£?<A<t<T}  >  2-5/2  (5) 

Now,  for  fixed  t  with  probability  greater  than  2-5/2  we  have 


B  A(t,&,z)  -  S  F~M(t,L,z) 


A <t .  .  <t  n  ..  h  .  . , 


-  In 


(3-  — - )'■«. 


By  which  implies 


h .  . 

i jr 


n . .  h... 


Hence,  the  left  hand  side  of  (6)  is  dominated  by 


9  *.9 
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•  ,  h  .  .  h  .  . 
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tji* 
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<7 

0  V 
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<7 

0  \ 
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nr 

(7) 


Now  there  exists  a  constant  A  such  that 
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J  -  X  <  (  »  +  ) 

tfc  n 


for  Vn . , >2. 
zt 


In  particular  choose 


,j  *  ( ' )  X  ql+2  (®  * ) 


We  then  have 


B  ASB  F-“ 


h.. 

<  l 


4  L  O 

v  r3 A<t .  •  <t  n  .. 

-  nr-  %tijr 


h .  . 


f~  Ai  [ *  ~~r~  dNnr8)> 

i  8 


(8) 


where  N..(s)  is  a  multivariate  counting  process  with  intensities 


N 


v.-fa)  -  l  exp  (a',  z  J  I(U  (z~ )  -  i)  (e)  JJt)  (9) 

*'«/  ymJ  *'V  V  V  •'J  V 


and 


M .  .ft)  -  N..(t)  -  J*  v  .  .fal  da 

1  7.  1  *1 


tj 


(10) 


are  orthogonal  integrable  martingales  by  Aalen's  (1975)  theorem  2.1 
on  page  14. 
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Let 


(11) 


X  .  .(a) 

tj 


!v  .  .  ( a  )  ,  n  .  >  2 

tj  *  i a 
0  .  n.  <  2 


Now,  consider  the  random  variable 


0-(w)  *  Ji  —s~ —  dN..(atu) 

1  l  n.  (<*>;  x° 
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Since  0.  is  positive,  we  have 


PI  In 

N  t0B  0j(ta)  *  0  if  we  can  show  that 


Ji;  e(02(^))  *  o 


E(0  (<*))  »  E  J 


ft  —J- -  dN  .  .( 8,vi) 

A  n.  (*>)  0 

18 


\*  .(a) 


L  ni* 

18 

ft  *  t 

because  M..  »  N..(t)  -}  \..(8)  da  and  }H(8,u>)  dM .  .(8,u)  are  square 

‘'J  Q  1'tl 

integrable  martigales  for  any  simple  process  H(8,u i)  as  defined  by 
Aalen  (1975)  on  page  8. 


We  have 


N  (v) 


M..(t)  *  S4j(t)  -  /  I  h..(a)  I(V  (a~ )=i)  J  (a)  vj,(e)  da 


A  v=2 


for  n.  >2 
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•  0  otherwise;  and 
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IJ  ■  M .7  V  *0  V 
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Plim 
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it 
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Since  0^t>0  by  A. 3.  Combining  equations  (8),  (13),  and  (16)  we  finally 
obtain  that 


Plim 

U-fco 


*  a  *  Z»  JL# 

BA(t,L,z)  -  S'  H(t,L,z) 


0  . 


(17) 


A  p  w 

In  other  words  our  estimate  B  is  asymptotically  equivalent  in  proba- 

A  M 

bility  to  Aalen’s  (1975)  estimate  BA(t, h,z).  But  the  consistency  of 
Aalen's  estimate  is  established  in  proposition  8.1  on  page  69  of  his 
thesis  (1975)  assuming  that  the  risk  sets  grow  large  like  N.  There¬ 
fore,  we  obtain 


Theorem 

*  F—M 

B  (t,b,z)  defined  by  equation  (2)  is  consistent  for  es¬ 
timating  B(t,b,z)  «  /*  v..(s,z)  de. 

A 

*F—M 

Weak  Convergence  and  Asymptotic  Variance  of  B  (t,A,a) 


From  the  martingale  convergence  theorem  given  on  page  29  and  the 
theorem  8.2  on  page  69  of  Aalen  (1975),  we  obtain  that 


✓s  [B*‘M(ttL,z)  -  B..(ttLtz) ) 

/  tu<»> 
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where  l^fe)  is  the  Weiner  process  and  /  is  a  normal 

process  with  independent  increments  where 
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Hence  the  asymptotic  variance  of  /JF  (ttA,z)  -  Bij(t,A,zJjj 
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Monte  Carlo  Comparlaon  of  Cox*s  Partial  vs.  Total  Likelihood 

The  performance  of  the  maxima  total  likelihood  estimate  end  the 
Cox  maximum  pert  lei  likelihood  estimate  are  compared  using  Monte  Carlo 
computer  simulation  techniques.  A  wide  variety  of  different  experimental 
conditions,  including  those  pertaining  to  an  actual  clinical  trial,  are 
used  to  evaluate  the  effects  of  changes  in  basic  parameters  on  the  else 
and  power  of  three  teats.  The  parameters  considered  Include  censoring 
(presence  or  absence),  the  baseline  hazard  rate,  the  number  of  covari¬ 
ates  and  their  distributions,  and  the  magnitude  of  the  covariate  effects. 
Distributional  data  for  the  estimates  of  the  covariate  effects,  includ¬ 
ing  confidence  bounds  for  selected  percentiles,  are  also  given. 

1.  Introduction 

In  recent  years  much  attention  has  been  given  to  the  problem  of 
accounting  for  covariate  information  in  survival  analysis  of  clinical 
trials  in  which  study  participants  may  have  been  censored.  Cox  (1972 ) 
addresses  this  situation  by  proposing  the  proportional  hazard  rate  model: 

Mt  |  s^P)  ■  XQ(t)  exp  [p  ,  z^]  , 

where  A(t  |  *i#p)  is  the  failure  rate  (hazard  function)  for  the  i-th 
patient,  is  the  covarlate  vector  for  that  patient,  and  p  is  the 
vector  representing  the  effect  of  the  covariates.  XQ(t)  is  an  arbitrary 
function  representing  the  baseline  failure  rate  —  the  failure  rate  for  a 
patient  with  covarlate  vector  -  0.  The  survival  function,  8(t  |  *t ,  p) 
and  the  failure  density  function  f (t  I  s^  ,  p)  are  then  given  by 

S(t  |  s^p)  mm{JQ  Y*lvp)d*)  "  (P*  **)  J  XQ(x)dx  j 

and 
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f(t|*1#p)  -  x(tUi#  p)  s(t|*i,  p) 

■  *Q(t)  exp  -|P'  ii  -  erp  (P'  z^)  J  XQ(x)dxj  , 

respectively. 

Much  analysis  has  been  done  in  the  exponential  case  where 
0 

XQ(t)  ■  e  0  (Glaaser,  6.967),  Sprott  and  Kalbfleisch  (1969)#  Breslov 

(1S?2)  and  Prentice  (1973))*  Simple  departures  from  the  constant 

hazard  model  are  accomidated  using  the  Weibull  model  where  XQ(t)  * 

0  +0  In  t 

e  (Prentice,  1973).  When  the  concern  is  chiefly  to  esti¬ 

mate  the  covariate  effects,  Cox's  partial  likelihood  (Cox,  1975)  is 
attractive,  for  it  ignores  the  formulation  of  ^(t).  Oakes  (1977) 
and  Efron  (1977)  establish  conditions  under  which  the  partial  likeli¬ 
hood  is  asymptotically  fully  efficient,  and  show  that  under  many 
circumstances  of  interest,  the  partial  likelihood  suffers  little  loss 
of  efficiency.  In  an  important  recent  work  in  this  area,  Lindsay 
(i960)  proposes  the  mixture  model  as  a  tool  for  evaluating  the  maxi¬ 
mum  partial  likelihood  estimates  in  the  covariate  cause,  and  gives 
conditions  under  which  it  is  fully  efficient  within  that  model,  as 
well  as  providing  a  lower  bound  for  the  asymptotic  variance  in  those 
cases.  Peace  and  Flora  (1978)  did  important  work  on  the  effectiveness 
of  various  tests  in  the  small  sample  situation,  using  Monte  Carlo 
methods  to  compere  the  size  and  power  of  tests  based  on  Cox's  maxi¬ 
mum  partial  likelihood  estimate  (MPLE),  $p ,  and  the  maximum  total 
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lilelihood  estimate  (MILE),  of  g ,  for  the  above  models  in  an 

experimental  situation.  In  this  paper  ve  expand  and  refine  their 
vork.  We  compare  the  performance  of  the  partial  and  total  likelihoods, 
using  three  test  statistics  based,  respectively,  upon  the  asymptotic 
normality  (Atf)  of  $ ,  the  likelihood  ratio  (LR)  criterion,  and 
the  asymptotic  normality  of  the  gradient  of  the  likelihood  function 
(DLOGL)  evaluated  at  the  null  hypothesis.  The  baseline  hazard  XQ(t)  is 
parametrised  using  zero,  one  and  two  of  the  nuisance  parameters  0  ,  respect. 

v 

ively.  The  test  statistics  are  compared  based  on  the  partial  likelihood, 
the  exponential  total  likelihood,  and  the  Weibull  total  likelihood. 


Our  results  are  divided  into  parts.  In  the  first  part,  ve 
consider  an  experiment  In  which  a  single  covariate  takes  the  values 
zero  and  one.  We  explore  the  size  and  power  characteristics  of  the 
LR,  AS  and  DLOGL  tests  for  the  MFLE,  the  exponential  MTLE,  and 
the  Weibull  MTLE,  and  how  they  change  due  to  censoring  or  variations 
in  the  null  hypothesis,  the  sample  size,  or  the  covariate  distribu¬ 
tion.  In  the  second  part,  ve  look  at  the  distribution  of  $  under 
these  circumstances.  In  part  three,  the  performance  of  the  various 
tests  is  evaluated  in  a  clinical  situation,  using  the  covariate 
distributions. and  censoring  patterns  found  in  a  clinical  data  base. 


2.  Test  Statistics 

« 

Let  denote  the  termination  times  for  the  N  test 
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«l 

participants  and  {x^}”  ,  ^  be  an  indicator  variable  vhere  x^  »  0 
If  the  participant  was  censored!  and  xi  *  1  if  the  termination  was 
a  failure.  The  log  of  the  total  likelihood  is  then  given  by: 

N 

log  1^,  »  ^  In  X0(t±)  +  In  ,  p)) 

L*1 

and  the  log  of  the  Cox  partial  likelihood  is  given  by: 

i»g  ic  -  2  *i(»'  zi  -  ^  I  3 ) 

i=l  JeR^) 

where  R(t^)  is  the  set  of  participants  at  risk  at  time  t^.  Using 

A  A 

a  modified  Newt  on-  Raphs  on  technique,  the  estimates  of  P  and  G  were 
found,  both  under  the  null  hypothesis  and  the  alternative.  Under  cer¬ 
tain  regularity  conditions  we  have,  as  N 


(l)  2 [log  L(P0)  -  log  L(p )]  ■  ■ — > 


where  k  is  the  number  of  parameters  being  tested  under  the  null 
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hypothesis.  Here,  p  is  the  vector  of  combined  parameters 

P  *  (P»©)  “  (P^jPg , • • .,5^# • • • ,Pn,©^» • • • j©^) 

where  0  parameterizes  >^(t).  Bje  hat  on  o  refers  to  the  maxi¬ 
mum  likelihood  estimate  of  o  and  the  zero  subscript  refers  to  the 
maximum  likelihood  estimate  under  a  specified  null  hypothesis.  Thus, 
if  the  hypothesis  is: 

VPi“PiO'  *  **><*•  tor 

versus 

K  :  &i  f  PiQ  for  seme  i  ,  1  <  i  <  K 

then 

-  A  A  A  A  . 

P0  *  P10*  *  ",Pk0,P(k*-l)0'  *  *  *,Pn0'®10'  ’  * ’'®m0' 

where  P^^  )q*  ••  *»^n0,®10'  *  *  ’*®o0  are  the  maxljDU®  likelihood  esti¬ 
mates  of  these  parameters  under  HQ. 

3.  The  Simulation 

Each  experimental  run  consisted  of  the  generation  and  analysis 
of  1000  samples.  The  survival  times  were  generated  randomly, 
according  to  the  formula 

-1  _ez< 

t4  -  S  (re  x  |  z  »  0  ,  p  -  0) 
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where  s’1{*  |  z  *  0 ,  p  ■  0)  is  the  inverse  of  the  survival  function, 
taken  when  p  ■  0  and  z  ■  0.  The  quantity  p  is  the  true  value  of 
P  far  the  simulation  and  S~^  is  generated  tabular ly  with  r  taken 
from  the  uniform  (0,  l)  random  number  generator; 


where 


«i  *  (n^  *  630  ,  360 , 016)  mod  (251  -  l)  . 


As  for  the  data  from  an  active  clinical  trial,  data  was  taken 
from  the  Western  and  Southeastern  Cancer  Study  Groups'  breast  cancer 
studies.  Three  covariates  were  chosen: 

-  treatment  type  (65 *3$  zeros,  3^-7*  ones)  , 

Zg  -  years  (y)  free  free  cancer  prior  to  entering  study  , 


•  0 

if 

0  <  y  <  1 

(U5.25&) 

«  1 

if 

l  <  y  <  5 

0*5-2*) 

-  2 

if 

5  <  y 

(  9.6*) 

«  0 

if  liver  metastases 

absent  at  entry 

(57.6*)  , 

-  1 

if  liver  metastases  present  at  entry 

Two  additional  covariates  were  defined  to  measure  interaction 
effects : 

Z5  -  (1  -  ZJ)(1  -  Z^  and  Zk  ■  (l  - 
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which  v«re  present  in  the  following  proportions: 


*4 

Percent 

0 

0 

0 

18.1 

0 

1 

0 

10.5 

1 

0 

0 

24.3 

1 

0 

1 

47.1 

The  partial  likelihood  estimate  of  0  over  the  breast  cancer 
study  ii  (0.47109  ,  -0.34805  ,  -0.26766  ,  -0.41*233)  when 

(z^ZgjZyZ^)  is  the  covariate  vector,  ©g  *  (-6.1771)  is  the 
nuisance  parameter  estimate  tinder  the  exponential  model,  and  9*  * 
(-6.9516 , 0.13988)  is  the  estimate  under  the  Weibull  model. 

The  method  of  censoring  is  also  derived  from  the  data  base. 
The  Kaplan-Meier  (K-M),  or  product  limit  estimate  of  the  censoring 
distribution  for  the  data  was  computed  and  stored  in  tabular  form. 

It  is  displayed  in  Figure  1.  A  random  censoring  time  was  evaluated 
for  each  study  participant  according  to  this  distribution  and,  if  it 
was  less  than  the  survival  time  for  that  participant,  the  survival 
time  was  replaced  by  the  censoring  time  and  the  participant  was  con¬ 
sidered  to  have  been  censored.  This  is  a  time  dependent,  and 
therefore,  covariate  dependent  form  of  censoring.  The  effects  of 
this  type  of  censoring  on  different  covariate  combinations  are 
discussed  later  in  the  paper. 


3.1.  The  Result*  -  Part  1 


Tables  1  and  2  shov  the  results  of  a  set  of  ninety  simulation 
runs  for  a  model  with  one  covariate,  z^  e  (0  ,  1} .  In  the  exponen¬ 
tial  case  0  ■  (0.0)  for  rows  1-6  of  each  table  and  0  •  for 
rows  7-9*  In  the  Weibull  case,  0  *  (1.0579*  1.0)  in  rows  1-6 
and  0  -  0*  in  rows  7-9*  In  the  Appendix  herein  we  show  that  the 
value  of  0  has,  in  the  absence  of  censoring,  no  bearing  on  the  value 
of  0  in  either  the  exponential  or  the  Weibull  models.  In  lines  1-7, 
z^  is  chosen  without  replacement  from  a  population  of  size  N,  con¬ 
taining  half  zeros  and  half  ones.  In  line  6,  z^  is  chosen  with 
replacement  from  the  same  population.  In  line  9,  is  chosen  with 
replacement  from  its  empirical  distribution  in  the  clinical  trial  data 
base. 

In  the  first  three  lines  of  each  table,  we  see  the  effects  on 
the  test  statistics  (under  :  0  ■  0)  of  variations  in  the  sample 
size  in  the  absence  of  censoring.  Power  curves  corresponding  to  the 
DLOGL  test  are  chosen  for  visual  comparison  and  are  displayed  in 
Figures  2  and  3>  Upon  inspection  of  either  the  tables  or  the  graphs, 
it  becomes  immediately  evident  that  the  partial  likelihood  is  much 
less  powerful  than  either  of  the  total  likelihoods  for  N  »  10.  The 
AN  test  performs  very  poorly  under  the  partial  likelihood,  the  LR 
test  does  better,  and  the  DLOGL  test  is  best  of  all,  though  it  still 
falls  well  short  of  the  same  test  in  the  total  likelihood  analysis. 

On  the  other  hand,  the  tests  in  the  total  likelihood  analysis  are 
very  similar  to  each  other,  although  the  DLOGL  test  still  appears 


ko- 


to  be  somewhat  more  powerful. 


Comparing  the  sizes  ot  the  tests,  we  observe  those  based  on 
the  Weibull  model  are  very  biased.  Using  the  Kolmogorov- Smirnov 
(K-S)  confidence  intervals  we  reject  the  hypothesis  with  95$ 
confidence  that  any  of  them  are  truly  A  r*Pre8entative  plot, 

using  the  DLOGL  test  is  given  in  Figure  7>  Under  the  partial  like¬ 
lihood,  the  DLOGL  test  is  again  not  j,  though  the  deviation  is 
not  so  great,  as  is  seen  in  Figures  4  and  6.  The  L.R.  test  appears 
to  be  "^2.)  ^  the  AN  test  is  borderline;  failing  the  K-S  test 
and  passing  it  once.  The  tests  based  on  the  exponential  model  appear 
to  be 

When  N  *  20,  the  partial  likelihood  compares  better  against 
both  of  the  total  likelihoods.  The  various  tests  based  on  the  partial** 
likelihood  are  much  closer  to  each  other  in  power,  and  much  closer  to 
the  corresponding  test  in  either  of  the  total  likelihoods.  The  sizes 
of  the  various  tests  are  stabilizing  as  well,  although  the  tests 
based  on  the  Weibull  model  are  still  rather  far  from 

By  the  time  N  ■  40,  the  irregularities  are  all  smoothed  exit, 
almost  down  to  random  fluctuations.  The  tests  are  now  almost  close 
to  w>d  the  powers  differ  only  slightly  between  tests  under  a 

given  model.  The  partial  likelihood  is  now  competing  wej_l  against 
both  the  exponential  and  Weibull  total  likelihoods.  The  DLOGL  test 
still  appears  to  be  the  strongest  test  in  all  caaes,  and  while  the 
AN  test  appears  to  be  slightly  stronger  than  the  L.R.  test  in  the 
total  likelihoods,  the  two  tests  are  almost  indistinguishable  in  the 


partial  likelihood,  with  perhaps  a  slight  preference  going  to  the 
L»S*  test* 

Lines  U  -  6  in  Tables  1  and  2  give  results  for  experiments 
identical  to  those  above,  except  that  they  correspond  to  testing 
Hq  :  P^  »  1  against  t-P^  f  1.  An  immediate  observation  is  that 
for  N  •  10  the  L.R.  and  AN  tests  have  no  power  at  all  against 
P^  >  1.  This  is  due  to  the  clustering  of  the  covariates,  so  that 
all  of  the  members  with  ■  1  fail  before  any  of  the  members  with 
-  0.  Under  these  circumstances,  the  MFX£  of  f)  is  infinite 
unler  the  alternative  hypothesis  and  the  analysis  falls  to  converge. 
Our  program  had.  no.  means  of  dealing  with  these  cases,  except  to 
throw  them  out,  which,  of  course,  seriously  weakened  these  tests  when 
it  happened  frequently;  namely,  when  P  was  large  and  N  was  small. 
The  DLOGL  test  computes  the  MFLE  only  under  the  null  hypothesis, 
therefore,  it  was  unaffected  by  this  problem.  Even  so,  it  too  was 
asymetrlc  in  P^  around  1,  being  more  powerful  against  P^  <  1 
and  less  powerful  against  P^  >  1.  This  bias  is  characteristic  of 
the  tests  based  an  the  partial  likelihood,  being  evident  for  N  -  20, 
40  as  well,  with  the  L.R.  test  being  least  affected  by  it.  On  the 
other  hand,  the  tests  based  on  the  exponential  total  likelihood  are 
not  visibly  affected  by  the  change  in  H^,  apart  from  random  fluctu¬ 
ations,  _ Just  as  they  did  before.  Tests  based  on  the  Weibull  total 
likelihood  showed  a  moderate  asymetry,  again  with  the  L.R.  test 
being  relatively  less  affected  than  the  others.  The  partial  likeli¬ 
hood  is  nearly  as  strong  as  the  Weibull  total  likelihood  on  the  left. 


though  significantly  less  strong  on  the  right,  while  the  Weibull  is 
weaker  than  the  exponential  at  all  points. 

Line  7  in  Tables  1  and  2  is  identical  to  line  3»  except  that 
censoring  is  imposed  using  the  K-M  estimate  of  the  censoring  dis¬ 
tribution  of  the  clinical  trial,  as  described  in  the  previous  section. 
When  z1  ■  0,  the  censoring  probability  is  .128  for  all  values  of 
P  ;  the  censoring  probability  for  *  1  varies  with  P  as  shown 
in  Table  3*  The  tests  all  lose  power  due  to  the  censoring,  more  so 
where  p  <  0 ,  due  to  the  greater  probability  of  a  sample  member 
being  censored,  than  where  5  >  0. 

Within  each  model  the  tests  retain  the  same  relative  strengths. 
Tests  based  upon  the  Weibull  MILE  suffered  relatively  less  loss  of 
power  than  the  tests  based  on  partial  and  exponential  total  likeli¬ 
hoods.  This  is,  perhaps,  due  to  the  effect  of  time-varying  hazard 
rate  on  the  censoring  distribution.  Examining  line  8  where  the  co¬ 
variates  are  selected  with  replacement  and  line  3  where  they  are 
selected  without  replacement,  we  observe  essentially  no  effect  on 
either  the  power  or  size  calculations.  However,  the  additional 
change  of  the  covariate  distribution  which  was  affected  in  line  9 
has  greater  influence.  The  various  tests  are  weakened,  but  not  uni¬ 
formly  so.  The  L.H.  test  is  weakened  more  on  the  right  than  on 
the  left,  while  the  AN  and  DLOGL  tests  are  weakened  more  on  the 
left  than  on  the  right.  The  bias  of  the  L.R.  test  is  observable 
on  all  of  the  models,  but  the  bias  of  the  AH  and  DLOGL  tests  is 
evident  only  on  the  MFLE  and  the  Weibull  MTLE. 


3.2  Remits  -  Bart  2 


Tablet  4  and  3  give  scaled  confidence  intervals  for  selected 
percentiles  of  $  when  N  ■  Uo,  under  conditions  similar  to  those 
mentioned  in  Part  1.  All  of  the  values  have  been  standardized  by 
subtracting  the  median  and  dividing  by  a  scaling  factor  equal  to 
the  Interquartile  range,  divided  by  1.3^9  (the  interquartile  range 

2  v 

for  the  normal  distribution  vlth  a  *  1). 

Deviations  from  normality  occur  sporadically  throughout  and 
are  attributable,  in  many  cases,  to  chance  fluctuations.  However, 
when  ?  ■  (l.O),  we  see  a  more  pronounced  trend  away  from  normality 
among  the  MFLE  and  the  Weibull  MILE.  The  right  tails  of  the  dis¬ 
tributions  are  generally  heavier  and,  in  the  Weibull  MTLE,  the 
median  is  biased  high.  However,  these  changes  are  not  discernable 
in  the  exponential  MTLE. 

3.3  Results  -  Bart  3 

Table  6  gives  the  results  of  tests  patterned  after  the  data 
taken  from  the  breast  cancer  data  base  mentioned  previously.  Samples 
of  size  40,  60  and  100,  were  taken  with  replacement  from  the  data 

base  itself  and  analyzed  both  for  size  at  P  *  0  and  power  at  P  »  (5* 
(restricted  to  the  number  of  covariates  being  analyzed)  against  the 
null  hypothesis  that  p  •  0.  The  covariate  vectors  used  are  (z^), 

(z^  ,  Zg  )  and  (z^  ,  Zg  ,  z  j  ,  z^ ),  with  the  distribution  of  the  terms 
given  in  Section  3.  In  the  exponential  case  •  ©£  and  in  the 
Weibull  case  5  m  0*  . 


The  performance  of  the  partial  likelihood  in  the  one-cow iate 
case  is  nearly  indistinguishable  from  the  performance  of  the  total 
likelihoods,  either  with  censoring  or  without.  The  DLOGL  test  is 
the  most  powerful  test  throughout,  while  retaining  good  size  charac¬ 
teristics,  and  the  L.R.  test  is  the  least  powerful.  Censoring  has 
a  minor  impact  on  the  power  of  the  tests,  but  in  this  situation  has 
no  effect  on  their  relative  power.  Even  at  N  *  100,  the  tests  are 
clearly  not  powerful  enough  to  decide  against  HQ  :  p  *  0;  larger 
samples  are  obviously  required. 

With  two  covariates,  the  tests  are  somewhat  more  powerful, 
but  unchanged  in  relationship  to  each  other.  Censoring  has  less  ef¬ 
fect  in  this  case  because  <  0,  which  reduces  the  probability  of 
censoring  for  many  sample  members. 

The  four-cowariate  case  is  not  strictly  comparable  with  the 
other  cases  due  to  the  interdependence  of  the  covariates.  This  can 
he  seen  in  the  power  of  the  vari*ms  tests,  which  is  lower  than  in  the 
two-oovariate  case  because  the  combined  effects  of  and  z^ 
nearly  cancel  each  other,  occurring  of  the  time. 

The  size  of  the  DLOGL  test,  and,  to  a  lesser  extent  the  AK 
test,  are  severely  affected,  especially  when  N  *  40.  Whether  this 
la  attributable  to  the  interdependence  of  the  covariate  valuea,  or 
whether  it  is  simply  due  to  the  large  number  of  covariates  and  their 
distribution  tsksa  independently  of  each  other,  remains  to  be  tested, 
although  the  latter  appears  likely. 

-1*5- 


k.  Computational  Note 


The  program  which  generated  these  results  was  written  in 
Fortran  IV  and  run  on  a  FDP  ll/5^  with  96k  of  core.  Bun  time  for 
each  simulation  varied  from  about  ^5  minutes  to  over  20  hours 
of  CHJ  time,  depending  on  the  complexity  of  the  analysis.  Well 
over  a  thousand  CHJ  hours  were  required  from  program  development  and 
simulation  run  time.  The  authors  wish  to  thank  the  Department  of 
Obstetrics  and  Gynecology  at  Cedars-Sinai  Medical  Center,  Los 
Angeles,  California  for  the  free  use  of  their  computer  for  this 
enormous  usage  of  computer  time.  This  work  would  have  been  impossible 
without  their  generous  aid. 

5«  Appendix 

Here,  we  show  that  when  the  co variates  are  time-independent 
and  censoring  is  absent,  then  the  MTLE  of  P  is  independent  of  0 
in  the  exponential  and  Weibull  models. 

Consider  first  the  Weibull  model, 

>>0(t|©0,«1)  «  exp  (©Q  +  In  t )  . 

Reparametrize  as  follows: 

/  ®i 

*  "  exp  ^  ©pn  )  «d  b  •  ®i + 1 

to  get 

-U6-  • 


\)(tla,b)  -  abt*’1  . 

Now,  suppose  t  arises  from  the  distribution: 

b. 

F(tllal*bi)  *  1  -  e*P-  la^t^  exp  (p*  .)]  . 

b 

If  (a^b^)  (eg  » bg ),  then  t^  -*  tg  such  that  a^t^  - 
b2 

■zV*  Now' 

N  ^  ml  b 

L(tllal*bl)  *  H  fai')i#ti1  exp[p'  *L  -  a^t^  exp  (p1  *  )1J  , 

1-1 

where  t^  and  z±  belong  to  the  i-th  member  of  a  study  of  size 
N. 


d  In  L(t^|a^  »  bx) 


^  tin  (a^b^)  -  (bA-l)  In  t^  exp  (p’  z^] 


*7  - 

J  ±ml 


m  L  [lij  -  *«  *i  hi  '*>  (6'  h>  • 


i-1 


Since  the  likelihood  is  maximized 


when  *1  “  b^  -  ^  f 


and 


i  *«  '  l  *«  H  Hi  (Si  H> 

i-1  i-1 


Substituting  in 


/  ®9  *>2 

‘a  '  (qViJ 


gives 


f  f1  i\  i/bii  , 

Z  *«  "  L  *l  \  )  hi  ejt*  V 

i-1  i-1  L  1 


We  now  have  the  following  identities; 


*2 

% 

% 


-  /  *2  1 


$  b 

X? 

bl 


The  results  for  the  exponential  case  nay  be  seen  as  a  special 
case  of  the  above  proof,  where  b1  •  -  bg  “^2  -1* 
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Legend  for  Tables 

Pq  :  value  of  p  under  null  hypothesis 
P  true  value  of  P 

Cov  :  method  of  selecting  covariate  values,  i.e., 

F  -  without  replacement  from  a  group  of  size  N 
containing  half  zeros  and  half  ones 
R  -  same  as  F  except  with  replacement 
E  -  with  replacement  from  the  values  of  z^  taken 
from  the  data  base 

Cen  whether  censoring  is  imposed  -  Y  ■  yes  and  N  *  no 
N  :  sample  size 

Test  :  L  -  likelihood  ratio  (  Lit  )  test 

A  -  asymptotic  normality  (Ah)  of  £  test 
D  -  asymptotic  normality  of  gradient  (DLOGL)  test 
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Size  and  Power  Results  for  Tests  Based  on 
the  Partial  and  Exponential  Cotal  Likelihoods 
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Distributional  Data  for  Partial  and 
Weibull  Total  Likelihoods 
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Table  6 

Size  and  Fewer  Assessnents  for  Test*  Based 
on  the  Partial,  Exponential  Total,  and  Welbull 
Total  Likelihoods,  Under  Experimental  Conditions 
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